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Abstract 



The ability of grains to rotate can play a crucial role 
on the collective behavior of granular media. It has 
been observed in computer simulations that impos- 
ing a torque at the contacts modifies the force chains, 
making support chains less important. In this work 
we investigate the effect of a gradual hindering of 
the grains rotations on the so-called critical state of 
soil mechanics. The critical state is an asymptotic 
state independent of the initial solid fraction where 
deformations occur at a constant shear strength and 
compactness. We quantify the difficulty to rotate by 
a friction coefficient at the level of particles, acting 
like a threshold. We explore the effect of this particle- 
level friction coefficient on the critical state by means 
of molecular dynamics simulations of a simple shear 
test on a poly-disperse sphere packing. We found 
that the larger the difficulty to rotate, the larger 
the final shear strength of the sample. Other micro- 
mechanical variables, like the structural anisotropy 
and the distribution of forces, are also influenced by 
the threshold. These results reveal the key role of 
rotations on the critical behavior of soils and suggest 
the inclusion of rotational variables into their consti- 
tutive equations. 



1 Introduction 

Granular media is the second most used material by 
mankind, only after water [TJ[2]. The annual produc- 
tion reaches almost 10 billions metric-ton/year. The 
industrial processing of these raw materials consumes 
nearly 10% of the total electric energy on Earth. 
Therefore, any optimization on the processing and 
any improve on the knowledge of GM will have a di- 
rect impact both on economy and welfare. Despite its 
ubiquitous nature, granular materials are far from be- 
ing fully understood, mainly because of the numerous 
and different behaviors that can be generated collec- 
tively from the interactions of the particular grains. 
Granular media can behave as a liquid, solid or gas, 
and its particular behavior depends on the way the 
system is stimulated. 

Soils represent a particular and important exam- 
ple of granular materials. They are the main subject 
of study for civil engineers dealing with settlements 
of buildings, earth pressure against retaining walls, 
stability of slopes and embankments, etc. Despite 
the rich complexity found in soils, like ratcheting and 
creep, there is a particular state that is independent 
of soil consolidation history, initial density and sam- 
ple preparation: the critical state [3J, SI El El EI IB] - In 
this steady state the deformation occurs at a constant 
average volume, void ratio, and effective stresses. 
The critical state has been characterized and stud- 



ied from long time; however, its microscopic origins 
remain elusive. 

The critical state is a fundamental concept in soil 
mechanics. It is linked with the so-called critical state 
constitutive models and failure criteria. Since the 
seminal work of Casagrande in 1936 [3], it has been 
found in sands and other types of granular media O 
IH G3 [91 Q33 ECQ. It is possible to reach the critical 
state through shearing, bi-axial or tri-axial tests. It is 
independent of the initial density, sample preparation 
or even shear rate, provided one is working in the 
quasi-static approximation. 

The values of the macroscopic variables charac- 
terizing the critical state depend both on external 
and internal variables. External control variables 
are the confining pressure and the shear velocity; in- 
ternal variables are the grains' size distribution and 
shape, the grain-to-grain friction coefficients, etc. It 
has been observed that the shape of the particles 
and its rotational freedom affect the final steady val- 
ues [H[S1[I1]. I n contrast with more complex shapes, 
spheres (or discs) favor rotations, because there is no 
geometric interlocking among them. If rotations are 
suppressed for spheres in some way, it is expected 
that the effective shear strength will increase. Suikcr 
and Flerk jS] showed by three-dimensional computer 
simulations with spheres that a complete hindering of 
rotations modifies the critical shear strength. Sim- 
ilarly, A. A. Peha, R. Garcfa-Rojo and H. J. Her- 
rmann [HE] studied the effect of particle shape on the 
critical state by using two particle kinds in 2D: almost 
isotropic and elongated. They found that both the 
void ratio and the coordination number at the crit- 
ical state are different (among other variables) for 
isotropic and elongated particles. Nevertheless, no 
systematic study has been done regarding the role of 
rotations as a continuous parameter, for instance, by 
hindering gradually the rotations. 

The main goal of the present work is to further 
explore the effect of restricting the rotation of the 
spheres on the macroscopic variables characterizing 
the critical state. This restriction is set by imposing 
a rotation threshold for each grain. The threshold 
is just the sum of the modulus of the normal forces 
on the grain times a rotational coefficient, which is 
the same for all grains. This procedure allows us 



to mimic the effect of geometric interlocking, but 
increasing the hindering to rotate in a continuous 
way. In other words, we are simulating the effect 
of complex shapes by adding a rotation restriction 
on rounded grains. The effective shear strength, the 
mean coordination number, the effective stresses, and 
the force and torque distributions, will be investi- 
gated as a function of this rotational threshold. The 
study will be performed through molecular dynam- 
ics simulations of a bi-dimensional and periodic poly- 
disperse system under simple shear. Section [2] intro- 
duces the critical state and the computational setup. 
Section [3] describes the discrete- element model, that 
is the interactions among particles and the integra- 
tion method for the equations of movement. Section[4] 
shows the main results and, finally, section [5] states 
the main conclusions from the work. 

2 The Critical State 

The critical state is defined in soil mechanics as 
the steady state where deformations occurs at con- 
stant average void ratio and effective stresses [131 [H] 

(Fig-[D- 

The definition of the critical state can be expressed 
mathematically as 

0(e) _ 9<p) _ %) 

at at dt ~ ' [ ' 

where (e) is the average void ratio e=Vy /Vs (i.e. 
the ratio between voids volume Vy and solids vol- 
ume Vg), (p) is the average isotropic pressure, (q) 
is the average deviatoric pressure and e q is the de- 
viatoric strain. In 2D, the mean pressure p is com- 
puted as p={u\ + 0-3) /2, and the deviatoric stress as 
q={a\ — (T3) /2, where a\ ((73) is the principal value of 
the stress tensor on the axial (tangential) direction. 

The inertia number / allows to quantify the degree 
of quasi-staticity of a given test [TH [TH [T71 HH] , and 
is defined as 

7=<W-^, (2) 

V Vwall 

where 7 is the shear rate, (d) is the mean diameter, 
p is the density of the solid grains and a wa u is the 
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Figure 1: (Color online) (Top) Sketch of the typi- 
cal responses for q, and (Bottom) of the volumetric 
strain, as functions of the deviatoric strain. 



Reference system 

Figure 2: Relevant variables for the interaction 
among two spherical grains. 



rection, rriij is the reduced mass, uj 
relative velocity, and 



Vj is the 



(4) 



confining pressure. In a quasi-static test, / < 1. In 
this work, I ~ 10~ 3 . 

3 Discrete Element Method 

We employed the Discrete Element Method (DEM), 
also known as Soft Particle method or Molecular Dy- 
namics, to simulate the direct shear tests on two- 
dimensional poly-disperse samples of spheres. Since 
the movements are 2D, both the translational and 
rotational degrees of freedom are integrated with the 
Velocity Verlet method Q1 CM [2TJ 122] • 

The normal force between two grains (see Figure [2]) 
is just the Hertz elastic law together with a damping 
term. It is given by 

f — i el 4- f da — h 
.In — J n ~Jn i 



-jmijivijofiij )h i j 
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where kij is the effective spring hardness, Rij is the 
reduced radius, hij = Ri + Rj — \ fi — fj | is the mutual 
penetration, 7 is the damping constant in normal di- 



For the tangential force we use the generalized 
Cundall-Strack model proposed by Luding [23], with 
a small correction for the tangential damping force. 
The algorithm is analogous for both sliding and 
rolling friction, and it reads as follows: As soon as a 
new contact appears (disappears) between two parti- 
cles, a virtual spring is created (destroyed) for each 
of the friction processes aforementioned. This spring 
stores the relative tangential displacement. Let us de- 
note the elongation of this tangential spring at time 
t by £ . The spring is always kept on the tangential 
direction by means of the transformation 



£ = £ - ("tj £ )fk. 



x £ x hi 



(5) 



By using the current tangential spring elongation, a 
temporary tangential friction force is computed 
as 

f t = -k t Rij( - rriijjtvt, (6) 

where kt is a tangential stiffness, and 74 is the tan- 
gential damping constant, which is introduced to im- 
prove the stability of the method [23] . Static friction 



is present whenever \$ \</J, s \fn\] otherwise we are 
on the kinetic friction regime. In the static regime 
ft=ft > an d the spring is updated as |* =£ 4- v t St. In 
the kinetic regime, the tangential force is given by 
ft=^d\fn\tij , where \x& is the kinetic friction coeffi- 
cient (typically fid = 0.8/z s ), and the tangential unit 
vector tij is defined dynamically as Uj — /j 3 / 1 fj* \ . 

The spring is updated as £,=^(ft + rriijlt'Vt)- The 
friction torque can be computed as 

r = kx f t , (7) 

where U is the branch vector from the center of the 
particle to the contact point. In total, four pa- 
rameters are required for each tangential interaction: 
&t>7t> A*s) but, in practice only [i s is needed, since 
the other parameters are chosen relatively to the nor- 
mal force parameters and fi s itself. 

A main goal of the simulation is to compute the 
macroscopic stress tensor <r a p . It is given by [24j ES] 

a aP = y X! fa l P J ( 8 ) 

cev 

where V is the representative volume where the av- 
erage is done, f a is a— component of the force, lp 
is the (3— component of the branch vector (the vec- 
tor joining the center of the particle to the point of 
contact) and the sum extends over all the contacts c 
among all particles inside the representative volume 
element. The structural anisotropy and mean coordi- 
nation number are extracted from the fabric tensor, 
defined as 

were n is the normal vector joining the centers of two 
contacting particles, and the sum extends over the 
contacts c. The structural anisotropy is a — 2(F 2 — 
Fi), where Fi is the i— th eigenvalue of the fabric 
tensor. 

Finally, we define a rotational threshold coefficient 
fi c to simulate the rotation resistance of complex 
shapes as follows: If the total torque \Ti\ on the par- 
ticle i is less than /j, c Ri ^ |^n,i|> where ^2 \Fn,i \ is the 
sum of the normal forces on particle i, then the total 



torque and the angular velocity are set to zero; oth- 
erwise they keep their current values. By this way 
it is possible to gradually hinder the rotation of the 
grains by changing the parameter \i c . Since the sum 
of normal forces magnitudes depends on the number 
in contacts, this mechanism could reflect the inter- 
locking of more general shapes, if rotation resistance 
is associated with geometric interlocking. 

4 Results and Discussions 

We have performed 2D simulations of periodic simple 
shear tests on a poly-disperse sample of 1024 discs. 
The radius of each particle is computed as |2Uj 

r= ~r — tttr — zt~r — v ( w > 

where z € [0, 1] is a random number. This particu- 
lar form is uniform in mass (no size dominates the 
mass of the system), and generates more small parti- 
cles than larger ones, as observed in natural soils. In 
contrast, the typical procedure of choosing the radius 
uniformly as 

R -^min ~l~ ^(^max -^min)) (H) 

will generate almost the same number of small and 
large particles and concentrate the mass of system 
on the larger ones. In the present simulations, the 
ratio i? max : ^?min = 4:1. The inter-particle static 
friction is 0.4. The rolling friction is 0.1. The full 
sample has been sheared (Figure |3b, assuring that no 
localization band has appeared. In the following, nu- 
merical quantities like 1.23(4) should be interpreted 
as 1.23(4) = 1.23 ± 0.04, i.e. the error is in the last 
digit, where error bars correspond to three times the 
standard deviation. We checked that /i c ~ 2 hinders 
the rotations almost completely, therefore there is no 
need to use larger values for this parameter. 

4.1 Macroscopic internal friction 

The ratio q/p as a function of the shear deformation 
e q is shown in Figure [4j First, the ratio grows almost 
linearly, until it gets the steady state, where it keeps 
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Figure 3: (Color online) Mean normalized horizontal 
translation as a function of the normalized height, 
for several values of the friction threshold /i c . L x and 
L z are the horizontal and vertical size of the sample, 
respectively. 

approximately the same mean value. The slope at the 
growing part represents a transient sample hardness. 
This hardness increases with the rotational threshold 
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Figure 4: (Color online) q/p as a function of the shear 
deformation 7 shear = Ax/L z , where Ax is the hori- 
zontal translation of the top wall, for several values 
of the threshold. The inset shows a zoom of the ap- 
proximately linear growing regime. 

As we stated before, the threshold fi c represents a 
difficult for the particles to rotate, analogous to the 



0.3 r 




0.5 1 1.5 2 



Figure 5: (Color online) Internal friction /x* for sev- 
eral values of the threshold. 

geometrical interlocking occurring in more complex- 
shaped bodies. It is well known [TJ] that irreg- 
ular bodies has a large macroscopic internal fric- 
tion angle /1* = T/<7 con f, where r is the shear 
stress and a con { is the confining pressure, because 
of geometric interlocking. The shear strength /j,* 
of the sample increases with increasing fi c (Fig- 
ure [5]). This support the interpretation of fi c as 
quantifying the geometric interlocking of more com- 
plex shapes. The internal macroscopic friction /z* 
grows like a saturated exponential with the form 
/i* = 0.292(l)-0.043(l)e" 12 ( 1 )^. If fx s increases, the 
final value of n* (equal to 0.292(1) in this case) also 
increases. Therefore, by tuning of ^ s and the thresh- 
old /i c , it is possible to model samples of spheres with 
very high values of /z*. 



4.2 Distributions of Forces and 
Torques 

The normalized distributions of forces and torques 
are presented in Figures [6] and [7] for several values of 
the rotational coefficient /i c . Unexpectedly, it barely 
affects the general shape of the distribution of forces 
and torques in the sample. 

It has been shown that normal forces, tangential 
forces and even the torques distribute as a power law 
for ranges smaller than the mean and as an exponen- 
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Figure 6: (Color online) Distribution of tangential 
forces at each contact, for several thresholds. In- 
sets: Zoom on given ranges. (Top) Log-log, (Bottom) 
Linear-Log 



tial for values larger than the mean [T5l , 



P(/)oc /(</>) ' (12) 
[e^-f/W, />(/) . 



For all the values of /i c , the exponents a and /3 are al- 
most the same, a = —0.16(1), while for larger forces, 
j3 = 0.93(2). Therefore, the forces' and torques' 
magnitudes distributions are barely affected by the 
threshold, although the relative importance of each 
particular distribution is slightly different. 
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Figure 7: (Color online) Sliding torque distributions 
for several thresholds. Inset: zoom in the range of 
small forces 

4.3 Mean Coordination number and 
structural anisotropy 

The mean coordination number, on one hand, is 
almost constant for all the values of /i c , with a 
value of (z) =4.1(2). On the other hand, the struc- 
tural anisotropy starts with a small value and decays 
rapidly to almost a constant value for large thresholds 
(Fig.j). 




Figure 8: (Color online) Structural anisotropy as a 
function of the thresholds. 



Although the mean contact number is practically 
the same for all thresholds, we found that the larger 



6 



the threshold the smaller the structural anisotropy. 
This result agrees with [T2J, where discs' and pen- 
tagons' samples were compared in biaxial tests and 
it was found that discs' samples have higher struc- 
tural signed anisotropies and smaller strengths, be- 
cause of smaller force anisotropy. If the threshold 
fi c can be associated with shape, then larger val- 
ues of fi c should reflect more irregular shape, higher 
shear strength and smaller structural anisotropy, as 
found. The higher shear strength with smaller struc- 
tural anisotropy can be explained by means of the 
forces' anisotropies, that reflects how the forces mag- 
nitudes are distributed along the contact network. 
Future works could check that increasing \i c enhances 
the force anisotropies and thus the shear strength. 
Since the structural anisotropy decreases rapidly to 
a constant value, the parameter /i c is shown to cap- 
tures some but not all the effects of irregular shapes. 
For increasing values of fx c , the connectivity of the 
contact network is kept constant in average, while 
the quantities that are carried over this network are 
changing, since (z) is almost constant while the struc- 
tural anisotropy decreases for increasing \x c . 

5 Conclusions 

In this work we studied the effect of a gradual hin- 
dering of the rotations on the critical state of a two- 
dimensional dry soil of poly-disperse spheres with size 
span 4:1. For this purpose, we have simulated a sim- 
ple shear test for several values of a rotational thresh- 
old friction coefficient fj, c , acting as follows: the rota- 
tion of a sphere is allowed only if the magnitude of 
the total torque on the sphere is larger than the sum 
of the magnitudes of the normal forces at all contacts 
times the threshold /i c . So, if ^ c = all torques (and 
its respective effect on rotations) are allowed, while 
if fi c — > oo, no rotation occurs. Actually, we found 
that values of /i c ~ 2 are enough to hinder almost 
completely all rotations. 

In order to characterize the effect of fi c on the 
macroscopic response of the system, five quantities 
were tracked: the stress ratio q/p, the total inter- 
nal friction fi*, the distribution of forces and torques, 
the mean correlation number (z), and the structural 



anisotropy of the fabric tensor. We found that the fi- 
nal stress ratio q/p is slightly increased for larger val- 
ues of /i c ; a consistent increment of the sample hard- 
ness before the critical state is reached can be easily 
observed. The normalized distribution of forces' and 
torques' magnitudes are barely affected in terms of 
its characteristic exponents, although is relative im- 
portance seems to decrease for higher values of fi c . 
The mean correlation number is barely affected by 
the rotational threshold fj, c , which means that the 
connectivity of the contact network is not changed. 
In contrast, the macroscopic friction coefficient //* 
does increase with /i c , like a saturated exponential, 
while the structural anisotropy decreases. This last 
result agrees with previous hindering strategies at the 
level of contact (instead of at the level of particles as 
here) |12) . This results enhances the interpretation 
of fi c as characterizing the effect of complex shapes 
and its natural geometric interlocking. Further works 
could check the evolution of the force anisotropies in 
order to explain the increasing shear strength with 
decreasing structural anisotropies. We expect the 
force anisotropies to increase with fj, c . The thresh- 
old [i c captures some but not all the features of more 
complex shapes. 

Different values of the static friction coefficient be- 
tween grains produce different values for the global 
internal friction coefficient fj,*, and the value of the 
later can be further increased by means of the thresh- 
old [i c . This relationship can be exploited to sim- 
ulate materials composed of spheres and with very 
high strengths, but further calibration with models 
including irregular shapes is needed. 

This work is a contribution to a better understand- 
ing of the critical state and the role of the rotational 
degrees of freedom. The exact dependence for dif- 
ferent macroscopic parameters, for different macro- 
scopic conditions, the role of forces anisotropies, and 
its quantitative comparison with more complex shape 
materials are topics of future work. 
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